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ABSTRACT 

We present general properties of ionized hydrogen (HII) bubbles and their 
growth based on a state-of-the-art large-scale (100 Mpc/h) cosmological radia- 
tive transfer simulation. The simulation resolves all halos with atomic cooling 
at the relevant redshifts and simultaneously performs radiative transfer and dy- 
namical evolution of structure formation. Our major conclusions include: (1) 
for significant HII bubbles, the number distribution is peaked at a volume of 
~ 0.6 Mpc 3 /h 3 at all redshifts. But, at z < 10, one large, connected network 
of bubbles dominates the entire HII volume. (2) HII bubbles are highly non- 
spherical. (3) The HII regions are highly biased with respect to the underlying 
matter distribution with the bias decreasing with time. (4) The non-gaussianity 
of the HII region is small when the universe becomes 50% ionized. The non- 
gaussianity reaches its maximal near the end of the reionization epoch z ~ 6. 
But at all redshifts of interest there is a significant non-gaussianity in the HII 
field. (5) Population III galaxies may play a significant role in the reionization 
process. Small bubbles are initially largely produced by Pop III stars. At z > 10 
even the largest HII bubbles have a balanced ionizing photon contribution from 
Pop II and Pop III stars, while at z < 8 Pop II stars start to dominate the overall 
ionizing photon production for large bubbles, although Pop III stars continue to 
make a non-negligible contribution. (6) The relationship between halo number 
density and bubble size is complicated but a strong correlation is found between 
halo number density and bubble size for large bubbles. 

Subject headings: cosmology: theory — large-scale structure of universe — early 
universe — galaxies: intergalactic medium — methods: numerical — radiative trans- 
fer 



INTRODUCTION 



Understanding the reionization process is a key challenge in cosmology. Present ob- 
servations provide extremely useful but still limited information with respect to the epoch 
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of reionization. On the one hand, the absorption spectrum observations of high redshift 
quasars from the Sloan Digital Sky Survey (SP SS) strongly suggest that the final reioniza- 



tion episode comes to complet ion at z ~ 6 (e.g.. iBecker et al.ll200ll ; iFan et al.ll2002l ; iBarkana 



2002 



Cen k McDonaldl |2002j) . On the other h a nd, the Wilkinson Microwave Anisotropy 



Probe (WMAP) observations (IPage et al.l 120061 ; ISpergel et al.l 120061 ) infer that the inter- 
galactic medium is already significantly reionized by some very early time, possibly z ~ 15. 
In combination, it s uggests that reionization process may be quite comple x and perhaps 
non-monotonic (e.g., ICeiJlioOa lliaiman k Holder] [2OO3I : IWvithe k CerJl2007l ). 



In addition, some othe r observations and analyses have also provid ed useful constraints 
on the reionization process. iTheuns et al.l (120021 ) ; lHui k Haimanl (120031 ) have shown that the 
observed temperature of the Lya forest at redshift z = 2 — 4 requires that the cosmological 
reionization occurred no earlier than z = 9 — 10, although He II reionization somewhat 
complicates this constraint. Analysis based on the SDSS quasar Strom gren sphere size 
suggest that the neutral hydroge n fraction at z ~ 6.3 is a few tens of percent (jWyithe k Loeb 
2004 ; iMesinger k Haimanl 120041 ). whereas analyses based on the evolution of Lyman alpha 
emitters from z = 6.5 to 5.7 im ply that a partially neutral IGM of neutral fraction of ~ 0.25 
is consistent with observations (IMalhotra k Rhoadal2004l; IStern et al.ll2005l ; lHaiman k Cen 



20051 ; IMalhotra k Rhoadsll2006l ). But the recent observations of Lyman alpha emitters begin 
to support the complete reionization at z ~ 6 as the SDSS quasar observations imply (e.g., 
Kashikawa et al.lhood : Lidz et al.lfcoOTl ). 



While a fully self-consistent r eionization picture is far from being painted (see lBarkana k Loeb 



200ll ; lLoebll2006l ; IFan et al.ll2006l . for a review), some key breakthroughs may rest in the red 
shift range z = 7—15, where upcoming observations, including 21cm radio and CMB observa 
tions, may be able to provide useful constraints. To understand the reionization process sev 
eral elegant analytical and semi-analytic models have been developed and predictions made 



(IMiralda-Escude et al. 



2005; 



Furlanetto et al 



2000 



Barkana k Loeb 



2006; Zahnetal 



20021; 



2007l : ICohn k ChaneJl2~0~0~7l ). These methods provide 



Furlanetto et al. 2004; Furlanetto k Oh 



ways to economically explore the large parameter space and have provided very important 
insights with respect to HII bubble evolution as well as the evolution of global quantities 
such as the ionization fraction. These methods, however, need to make necessary simplify- 
ing assumptions, which would limit the scope of their applicability or the accuracy of the 
predictions. 

We will take a complementary approach by making detailed numerical simulations with 
fewer assumptions. In earlier simulations, a choice is forced to be made between simulating a 
large volume with very limited re s olution or a sufficient resolution with too sma l l volume (e.g., 
GnedirJl2~00~ol ; Iciardi et al.l I2OO0I : iRazoumov et~aD I2OO2I ; ISokasian et al1l2003l . l2004j ; iGnedin 
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2004 ). However, both a large volume (~ 100 Mpc/h) and an adequate resolution are neces- 
sary in order to follow the reionization sources and sinks properly. A lar ge simulation box 
is required because the sources in question are extremely highly biased (IBarkana fc Loeb 
20041 : iTrac fc Cenll2006h and the bub bles have sizes of tens of Mpc prior to complete per- 



(e-g., 



Furlanetto et al.l 120041 ) . whereas resolving the bulk of sources of halo mass 

10 6 M© and spatial resolution of a few 



eolation 

~ 1O 8 M dictates, at least, the mass resolution of 



kpc. With a uniqu e hybrid (dark matter + baryons + stars + radiation) computational code 
(ITrac fc Cenll2006l ) we have crossed the threshold to be able to simultaneously simulate a 
large volume and have an adequate resolution to identify the bulk of the ionizing sources at 
high redshift as well as to have an approximate yet ade quate treatment of radiation transfer . 
Recent direct simulations performed by other groups (llliev et al.l l2006at IZahn et al.l 120071 ) 
have yielded very important results but the inability to adequately resolve dark matter halos 
of ~ 10 s M renders the results uncertain. Additional features of our method include (1) 
following simultaneously the evolution of structure formati on and radiative t ransfer, instead 
of performing radiative tra nsfer as a post-proces sing; step (llliev et al.l l2006bl ) or adding un- 



resolved halos analytically (IMcQuinn et al. 



20071 ). (2) a self-consistent, albeit still uncertain, 



treatment of metal-enrichment process, hence the spatially varying transition from Popu- 
lation III (Pop III) to Population II (Pop II) initial mass function (IMF). Therefore, our 
simulations have more direct numerical treatments of radiation sources, small-scale clump- 
ing, and self-shielding than both post-processing and semi-analytic models. 

In this paper, we analyze the evolution of HII bubbles in a 100 Mpc/h simulation 
box. In particular, we study how HII bubbles grow and what physics determines the growth. 
Understanding the evolution of HII bubbles is closely related to predicting fu ture observation 



results of the redshifted 21cm line and cosmic microwave background (see iFan et al.l 12006 



for a review). Among large volum e reionization simula t ions the morpho l ogy o f HII regions 
were studied in llliev et al.l J2006ah : IZahn et all J2007h : IMcQuinn et aD (l2007h . In §2, the 
details of our simulation are described. Results of simulations are given in §3, followed by 
discussion and concl usions in §4. We use the following cosmological parameters based on the 
WMAP3 results (see lSpergel et al.ll2006l . and references therein): Q m = 0.26, = 0.74, fi& 
= 0.044, h = 0.72, cr 8 = 0.77, and n s = 0.95. Throughout the paper both length and volume 
are given in comoving units. 



2. SIMULATIONS 



Our simulation is generally based on the numerical methodology described in lTrac &: Cen 

( 120061 ) and similarly utilized an N-body algorithm for dark matter, a star formation prescrip- 
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tion, and a radiative transfer (RT) algorithm for ionizing radiation. However, we have taken 
some simpler approaches in this initial step to satisfy the computational challenge of large 
volume, high resolution simulations of cosmic reionization. In particular, we use an alterna- 
tive RT algorithm and do not use a halo model for prescribing baryons and star formation. 
Here we summarize the main components and describe the modifications. 

In a L = 100 Mpc/h simulation box, a high resolution N-body calculation for 2880 3 
dark matter parti cles on an effectiv e grid with 11520 3 cells is performed using a particle- 



multi-mesh code ( jTrac &: Peru 120061 ) . The particle mass resolution is 3.02 x 10 M© and 
approximately 33 particles make up a lO 8 M /h halo. The comoving grid spacing is 8.68 
kpc/h and approximately 12 cells make up the virial volume of a 10 8 M /h halo. We identify 
dark matter halos in post-processing rather than during the course of the simulation. 



We do not use the halo model of iTrac fc Cenl (120061 ) for prescribing baryons and star 



formation. Instead, an alternative approach is taken where we calculate the local matter 
density p m and velocity dispersion a v for each particle. The baryons are assumed to trace the 
dark matter distribution on all scales and we obtain the local baryon density p& = p m (Qf,/fl m ) 
and gas temperature T = pa^/(3k). Star formation is only allowed to occur in particles with 
densities p m > 100p cr i t (z) and temperatures T > 10 4 K. This cut in the temperature-density 
phase-space effectively restricts star formation to regions within the virial radius of halos 
where efficient atomic line cooling allows the gas to dissipate energy and further collapse 
to very high densities. We also differentiate between the first generation Pop III stars and 
the second generation Pop II stars by following t he chemical enrichm ent of the interstellar 



medium and intergalactic medium as described in ITrac fc Cenl (120061 ) . 



Our radiative transfer algorithm for ionizing radiation is based on a photon-advection 
scheme which is much less computationally expensive than ray tracing. Particles, sources, 
and photons are collected on a RT grid with 360 3 cells. Each cell spans 278 comoving kpc/h 
and the RT time step is set by the light-crossing time. For a source cell, the excess source 
photons are propagated to the 26 neighboring cells. The advection is photon-conserving and 
the isotropic redistribution uses a weighting function that is proportional to 1/r 2 . However, 
for a non-source cell with excess radiation, the advection is generally anisotropic. For a HII 
region, photons originating from a central source propagate in the direction coinciding with 
decreasing radiation flux. Therefore, we propagate photons by looking for gradients in the 
radiation field. Consider a non-source cell with cell indices (i,j,k) and radiation density 
n 7 (i, j, k). Radiation can propagate to an adjacent cell with indices (i + di, j + dj, k + dk) if 
either of the gradient terms, 

A + = n 7 (i + di,j + dj,k + dk) -n 7 (i,j,k), (1) 
A_ =n 7 (i,j,k) -n 7 (i-di,j -dj,k-dk), (2) 
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is negative. The first gradient term indicates the downstream or expansion direction while 
the second term indicates the upstream or source direction. For the 26 possible neighbor cells, 
we count how many cells satisfy the above criteria and redistribute the excess photons equally 
among them. If none of the 26 neighbor cells satisfy the above criteria, then the central cell 
is a convergent point and we redistribute the photons equally among the neighbors. 

For an isolated HII region with one central source, if one of gradient terms is false, the 
other is generally false too. However, this may not be the case near the interface of merging 
or overlapping HII regions. If a weaker radiation field is trying to expand into the vicinity 
of a stronger radiation field, then the first gradient term will be positive even though the 
second gradient term will be negative. Therefore, it is necessary to use either of the gradient 
terms to determine the direction of radiation propagation . In the app e ndix, we compare the 
photon-advection scheme with the ray-tracing scheme of iTrac fc Cenl ( 120061 ) and show that 
this simpler approach correctly captures the radiative transfer for a significantly majority of 
the reionization epoch. 

For each particle, we store 12 floating-point variables: three dimensional position, three 
dimensional velocity, matter density, baryon density, temperature, stellar fraction, HI frac- 
tion, Hel fraction, and Hell fraction. We calculate the ionization and recombination for each 
particle individually rather than using the lower resolution density field defined on the RT 
grid. This allows us to correctly account for the clumping factor and self-shielding of dense 
gas down to small scales of several comoving kpc/h. 

The simulation was run at the National Center for Supercomputing Applications (NCSA) 
on a shared-memory SGI Altix with Itanium 2 processors. We used 512 processors, 2 TB 
of memory, and approximately 80000 cpu hours. With nearly 24 billion particles, this is 
the largest cosmological N-body simulation run to date. This is also the first reionization 
simulation with a L = 100 Mpc/h simulation box that can directly resolve dark matter halos 
down to virial temperatures of 10 4 K. 

In pos t-processing, we ha ve identified dark matter halos using a friends-of-friends (FoF) 
algorithm ( Davis et al.l Il985l ) with a standard linking length of b = 0.2n -1 / 3 , where n is 
the mean particle number density. Figure 1 compares our mass functions with the ana- 



l ytic p rediction of iPress fc Schechterl (119741 ) and the empirical prediction of IWarren et al. 



(120061 ). For redshifts z < 10, our results are in very g ood agreement wi t h recent works on 
the mass function of high redshift dark matter halos (IReed et al.l 120071 ; iLukic et al.l 120071 : 
Cohn fc Whitd 120071 ) . For higher redshifts z > 15, we systematically under-resolve halos 
because our starting redsh ift of z = 60 is too late to capture the nonlinear gravitational 
collapse. IReed et al.l (120071 ) have suggested that simulations must start ~ 10 — 20 expansion 
factors before the redshift at which results converge. We have found that we can accurately 
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capture the format ion of halos at z 
dTrac fc Cenll2006h . 



15 — 20 with a smaller simulation starting at z — 300 



Our star formation pre scription resembles those used in hydrodynamic simulations (e.g., 
Springel fc Hernquistll2003l ) and we obtain very similar results. In Figure 2, we compare the 
star formation rate p(z) with an analytical model calibrated using hydrodynamic simulations 
( iHernquist Springelll2003l ). corrected for our cosmology. The overall shape is in very good 
agreement at all relevant redshifts. At z = 6, our amplitude is 1.25 times larger and this 
difference is simply due to the fact that we have chosen a coarse value of c* = 0.06 for 
the star formation efficiency For our purposes, the overall amplitude of the star formation 
rate is unimportant since it is degenerate with the radiation escape fraction. We have 
correspondingly chosen a radiation escape fraction of / csc = 0.15 in order to have the redshift 
of complete overlap occur at z ~ 6. 



3. RESULTS 

The first reionization region appears when regions of dense baryons turn on star forma- 
tion. Once the radiation sources begin ionization around them, the global ionization regions 
appear as shown in Figure 3. At the early stage of reionization, isolated ionization regions 
are easily found (see z = 11.228 in Figure 3). As more radiation sources develop, the isolated 
bubbles get connected together along filaments. The process of HII bubble merging is very 
complex and difficult to treat without detailed simulation. Visually, in Figure 3, the com- 
puted HII bubbles do not appear to be close to spheres, as will be confirmed quantitatively 
below. After z ~ 10 one large connected network of HII bubbles begins to dominate the 
simulation volume, the HII percolation process enters the "cannibalistic" phase where the 
dominant HII bubble rapidly swallows other HII bubbles, as evident in the z = 9.3, and 7.2 
panels. In the following subsections, we will present quantitative results. 



3.1. Size of HII bubbles 



3.1.1. Characteristic size 



The size of HII bubbles is a basic quantity that has received significant attention 



(Fur 



anetto & Ohll2005l : iFurlanetto et al.ll2006l : llliev et al.ll2006al : iMcQuinn et al.ll2007l : IZahn et al 



20071 ). We define HII bubbles as follows. First, we mark all the cells that are ionized at 
> 50% level. Then, a bubble is defined by all such cells that are connected at least by one 
face; this is practically done by grouping cells using a friends-of-friends (FoF) algorithm with 
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a linking length equal to simulation data ouput cell size, i.e. ~ 0.14 Mpc/h. In addition, at 
least two cells are needed to make up a single bubble. 

Since the shape of HII bubbles is not close to a sphere as shown in §3.2, we present, in 
Figure 4, the volume distribution of the found HII bubbles instead of the radius of a sphere. 
Note that in previous analytical studies HII bubbles are often assumed to be s pherical (e.g., 



Pritchard et al.l 120071 ) or have a characteristic size (e.g., iMcQuinn et al.l 120071 ). We present 



the size distribution of z = 5.7, 6.3, 8.9, 13.5, 17.2, and 20.6. 

In this paper, we define the characteris t ic size of bubbles based on their number fraction 



not volume fraction as used in llliev et al.l (j2006al ). Therefore, our characteristic sizes well 



describe th e exist ence of many small non-connected bubbles. Meanwhile, the definition of 



Zahn et al.l (120071 ) is affected by volume occupied by HII bubbles because it measures an 



ionization fraction within a certain smoothing radius. 

The volume distribution by the number fraction shows that there are three characteristic 
sizes of the HII bubbles before the simulation box is dominated by a single bubble. The 
characteristic volumes are 0.6, 0.03, and 0.006 Mpc 3 /h 3 . These characteristic volumes do 
not change even though more regions become ionized. But the least volume of the possible 
bubbles is ~ 0.005 Mpc 3 /h 3 because our bubble finding method defines the smallest bubble 
as at least two connected simulation cells. Therefore, the smallest-volume peak in the HII 
bubble distribution may be limited by our radiation cell size and should be considered as an 
upper limit for the peak at the smallest scales. 

Other two characteristic sizes of ~ 0.03 and ~ 0.6 Mpc 3 /h 3 are real features of the size 
distribution. The peak at ~ 0.03 Mpc 3 /h 3 represents HII bubbles that are not limited by 
our radiation transfer resolution and should be real. The peak at ~ 0.6 Mpc 3 /h 3 proba- 
bly represents typical mature HII bubbles produced by individual galaxies (probably with 
their satellite galaxies), which may subsequently merge with other bubbles. This peak is 
maintained at all times until when the percolation permeates the entire volume to complete 
the reionization, as seen from the difference between z = 6.3 and z = 5.7 in Figure 4. The 
properties of these characteristic sizes will be further analyzed and understood in §3.1.3 and 
§3.1.4. 

The appearance of few dominant large bubbles is easily found by the change of volume 
fraction for different sizes of bubbles. As shown in Figure 4, the volume fraction is dominated 
by few large bubbles after z = 13.5. This result is also visually recognized in Figure 3. In 
the plot of z ~ 11.2, large bubbles begin to appear while many small bubbles around newly 
formed sources dominate the number fraction in Figure 4. Therefore, the characteristic sizes 
of 0.6, 0.03, and 0.006 Mpc 3 /h 3 can be maintained until the scale of large bubbles is big 
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enough to percolate small bubbles that are distributed between large bubbles. 



The volume distribution of HII bubbles in llliev et al.l (j2006al ) can be compared to our 
results. Both simulations use the same definition of HII bubbles and the bubbles are found 
by FoF method. The first difference is that in our simulation large bubbles of > 10 4 Mpc 3 /h 3 
app ear later th a n z ~ 15 when the global volume-weighted ionization fraction is about 0.04. 
But llliev et al.l (]2006al ) found the existence of those large bubbles around z ~ 15 when in 



their simulation the global volume-weighted ionization fraction is 0.05. This difference is 
also related to different characteristic sizes that found in our simulation. In our simulation, 
most HII bubbles are much smaller than > 10 4 Mpc 3 /h 3 . Even the found characteristic size 
of intermediate-size bubbles is ~ 0.6 Mpc 3 /h 3 . llliev et al.l (j2006al ) generally shows larger 
bubbles for all ranges of the global ionization fraction and the formation of large bubbles 
at the earlier times than ou rs. Moreov e r, the largest fraction of bubbles is explained by 
intermediate-size bubbles in llliev et al.l ( ]2006al ) while our results show the dominance of 
small bubbles for almost all simulation time. 

We think this is pro bably due to the d ifference in the adopted cosmological model in 
which the model used by llliev et al.l (j2006al ) has a higher erg and the fact that they require 
the universe to be ionized much earlier such that their compute d r fi matche s the r esults 
from the first-year WMAP results. The low mass resolution of llliev et al.l (j2006al ) also 
does not permit the formation of sma ll bubbles because it cannot capture the formation of 
low-mass DM halos (IZahn et al.l 120071 ) . Finally, the implementation of baryon physics like 
ionizing sources in their simulation is quite different from our prescription given in §2. The 
importance of this difference is explained in §3.1.3 and §3.1.4. 

In Fig u re 4, we also present the bubble size distribution based on the definition of 



Zahn et al. 



(120071 ) for comparison. The same bubble boundary threshold of 90% is adopted 



i n this measurement. For similar global volume-weighted ionization fractions as IZahn et al. 



( 120071 ). the size distribution from our simulation has its peak at the smaller radius than their 
size distribution. The existence of many resolved small bubbles re sults in the differe nce of 
the size distribution even though we adopt the same definition as IZahn et al.l (120071 ). For 
example, we find two peaks of size distri bution at ~ 0.1 an d 0.5 Mpc/h for the volume- 
weighted ionization fraction of 0.12, while IZahn et al.l (120071 ) found a single peak at ~ 0.7 
Mpc/h for the volume-weighted ionization fraction of 0.11. Comparing the results for the 
volume-weighted ionization fraction of 0.52, our distributio n shows a single p eak size, but the 
size is ~ 2.5 Mpc/h that is smaller than ~ 4.5 Mpc/h of I Zahn et al.l (120071 ). In particular, 
the resolved small bubbles makes the peak probability of the size distribution, i.e. the height 
of the peak, insensitive to the change of a global ionization fraction. Therefore, even with 
the same definition of the bubble size, our result shows a slightly different result although a 
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general result is consistent with each other. 



3.1.2. Bubble merger history 

As it turns out, the percolation of bubbles is also related to the transition of major 
ionizing photon production from Pop III to Pop II stars. 

To shed light on the overall reionization process, we follow the change of ionizing envi- 
ronment at four different location in the simulation box that are occupied by four individual 
bubbles that all started at z ~ 21, shown in Figure 5. In the plot, we track the size change of 
bubbles that occupy the four locations. We see that these bubbles initially, formed at about 
the same time but at four different locations, display a variety of histories, clearly due to the 
rich, disparate structure formation histories and the evolution of ionizing environment. One 
striking and very clear signature is some intermittent, rapid downturns in the HII volume 
during an individual location's evolution history. This is due to a rapid transition of the 
IMF from Pop III to Pop II in the bubble, which in general occur at the redshift range of 
z = 10 — 15. This is consistent with the overall contributions of Pop II stars and and Pop III 
stars as a function of redshift shown in Figure 6, where we see that ionizing photons from 
Pop II stars begin to become more dominant between z = 10 and 15. But this transition 
occurs spatially at different time such that each bubble position has a slight different history. 
We note that this transition time is not the same as the epoch when star formation rate is 
dominated by Pop II population, since Pop III stars are more efficie nt producers of ionizing 



photons than Pop II stars for a given amount of star formation rate. iMcQuinn et al.l (120071 ) 
pointed out that the nature of the ionizing sources plays an important role in shaping the 
bubbles and our results quantitatively confirm their conclusion by being able to perform a 
detailed calculation that allows for a spatially resolved transition from Pop III to Pop II 
stars. 

Merger history of the bubbles is also consistent with the changing total number of 
bubbles and the dominance of the largest bubble over small bubbles after z < 10. As 
percolation of the bubbles proceeds, the total volume of HII bubbles increases while the 
number of individual bubbles decreases as shown in Figure 7. That is, percolation of the 
bubbles becomes so important that the largest network of connected bubbles dominates the 
simulation box. Although the number of small bubbles also begins to decrease after z < 10, 
the number fraction of small bubbles increases because the size of the small bubbles is not 
large enough to be connected to other bubbles (see Figure 4). 
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3.1.3. Dark matter halo and bubble size distribution 

Basically, the clustering of radiation sources affects the percolation domination and 
the spatially varying different transition epoch from Pop III star formation to Pop II star 
formation. Because the two small characteristic size peaks found in Figure 4 are maintained 
by the formation of new radiation sources, the correlation between the bubble sizes and the 
DM halos can be expected to explain a physical reason behind the bubble size distribution. 

In Figure 8 we show the number density of halos as a function of HII bubble volume. All 
DM halos whose centers are inside the bubble are said to belong to the bubble. We separate 
the member DM halos to three separate mass ranges, M ha i < 10 s M /h, 10 8 M /h < M ha i < I0 9 M /h, 
and 10 9 M Q /h < M ha i . We see that there is an upturn of halo number density for all mass 
ranges for the smallest bubbles. This indicates that these bubbles are rapidly expanding 
due to newly formed sources while the number of halos stays approximately constant. These 
smallest bubbles do not contain DM halos of Mhaio > 1O 9 M , corresponding to the peaks 
of ~ 0.11 and 0.19 Mpc/h as found in §3.1.1. The halo number density decreases, as one 
moves to the right, and then levels out at ~ —1 to of the displayed x-axis. The most 
interesting feature is that there is a sudden rise of halo number density of all masses toward 
the largest bubbles, and points to the direction that percolation is expected to continue 
with time, consistent with our analysis of the actual bubble size evolution. The boundary 
between the percolation dominance and the source formation dominance is the characteristic 
volume of ~ 0.6 Mpc 3 /h 3 . This feature is described as the beginning of the flat DM halo 
density distribution in Figure 8. A large drop in the minimum number density is caused by 
difference of bubble age among DM halos of the same mass. For example, the same mass 
DM halos of 10 8 ' 5 M /h collapse at different epoch and it results in the differently matured 
HII bubbles around them. So, the number density of DM halos can have a large dispersion, 
in particular, for the intermediate-mass DM halos. 

As characteristic sizes of HII bubbles do not change much until z ~ 6, the DM halo 
distribution does not change as shown for z ~ 8 and 10. So, we find that the size distribution 
of HII bubbles are mainly regulated by an overall mass density that determines the evolution 
of baryons. The latter property indicates highly biased galaxy formation and suggests that 
stellar mass density is expected to rapidly rise in high density regions. Whether this property 
may be preserved requires more detailed analysis. 
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3.1.4- Dependence of bubble size on ionizing photon production 

The rate of ionizing photon production per unit volume is likely the main physical factor 
to determine the evolution of HII bubbles, while the broad correlation between star formation 
rate and DM halo density has likely produced some of the interesting features presented in 
the previous sub-section, because the latter is essentially proportional to the total, integrated 
ionizing photon number density. An ionizing photon production rate inside HII bubbles is 
estimated to be proportional to star formation rate times the ionizing photon production 
efficiency that is a function of the IMF. Figure 9 shows the ionizing photon production rate, 
separately from Pop III stars and Pop II stars, as a function of the bubble size at z = 8 
and z = 10, respectively. The flat part of the average ionizing photon production rate 
represents the value which is approximately set by the average baryon number density and 
recombination rate inside the bubbles. It is clear that for the smallest bubbles with size 
< 0.3 Mpc/h Pop III stars dominate the ionizing photon production rate at both z — 10 
and z — 8, indicating that these bubbles are produced by the first generation of stars formed 
in those locations which are likely relatively removed from large halos where star formation 
has gone through more than one generation and the gas has been significantly enriched with 
metals, consistent with Figure 8. We expect that Pop III galaxies may be present in small 
bubbles until z ~ 6. For the large bubbles (size > 10 Mpc/h) the contributions from Pop 
III stars and Pop II stars seem balanced at z ~ 10. But, by z ~ 8, these large bubbles 
begin to be dominated by Pop II stars, although Pop III stars continue to make significant 
contributions of ionizing photons. 

We show the careful consideration of this transition is needed to understand the reason 
why we find the bubble size distribution given in Figure 4. We believe that the overall 
evolution of HII bubbles depends on detailed modeling of Pop II and Pop HI stars. A 



neglec t or non-detailed treatment (e.g., IZahn et al.l 120071 ; iFurlanetto fc Ohll2005l ; llliev et al. 



2006bj) of the Pop III stars would translate to a significant change in the characteristics of 



bubble evolution. 



3.2. Shape of HII bubbles 

The shape of HII bubbles are far from a sphere unlike what is assumed in most analytical 
studies. Our result is the first effort to present a quantitative measurement of HII bubble 
shape in the research of cosmic reionization. In order to quantify the shape of the bubbles 
as ellipsoid, first, we calculate the inertia tensor: 
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where x is the coordinate of every cell that is included in each HII bubble, and x c is a 
geometrical center of the bubble. From the above defined Iy, we obtain the square roots of 
the three eigenvalues, a > b > c. We measure the shape of the bubbles that have at least six 
member cells. Therefore, the smallest volume of the bubble is ~ 0.016 Mpc 3 /h 3 . The shape 
is expressed as two parameters e± and e 2 : 



b 2 / c 2 

e ' = V 1 -?' e2 = V 1 -r (4) 

After percolation becomes a main process in growth of HII bubbles, the above measurement 
of shape is useless. At late stage of reionization, the complicatedly connected bubbles are not 
similar to ellipsoids. Moreover, neutral hydrogen regions can be surrounded by HII bubbles 
such as a grape-like structure. Therefore, we derive the t\ and e 2 for bubbles before the 
largest bubble is comparable to the simulation volume. 

The bubble volume-weighted shape distribution is estimated in the following way: 

P( ei ,e 2 )= f P(e u e 2 \V)P(V)dV. (5) 



Jv 

Figure 10 shows the concentration of volume-weighted HII bubbles in (ei, e 2 ) plane. Note 
that a spherical bubble would have e% = e 2 = 0. In Figure 11 we show the bubble shape as 
a function of bubble size, with the largest bubble separately shown as the star symbol, for 
z = 13.8 and z = 10. The shape of the bubbles follows a complicated dependence on their 
sizes. The shapes of small bubbles are, however, less accurately calculated because of limited 
resolution. Clearly, most of the volume in HII bubbles is far from spherical. The fact that e\ 
is quite close to 1 and e 2 has a wide range indicates that bubbles tend to have shapes ranging 
from a cigar to a ruler. Note that at z < 12 the shape distribution is basically dominated 
by the largest bubble. For example, at z = 10.0, the most probable shape parameter is 
close to ei ~ 0.9 and e 2 ~ 0.6 of the largest bubble, corresponding to axial ratios b/a = 0.4 
and c/b = 0.8. This high complex non-sphericality would introduce substantial inaccuracies 
in calculations based on analytic modeling of observables of the high-redshift universe (see 



Fan et al.l 120061 . for a review). The complicated shape of HII bubbles is largely determined 
by the clustering of ionizing sources, as bubbles grow around large-scale structure, primarily 
filamentary networks. Figure 12 gives the three orthogonal projections of a single randomly 
chosen bubble which show characteristic of bubbles resulting from mergers of smaller bubbles. 



3.3. Non-Gaussianity of reionization 



We sample 5000 spheres randomly placed in the simulation box and measure the mass- 
weighted ionization fraction of the spheres. That is, we measure xtubbie = (mass of HII) / 
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(mass of H) for each samp led sphere. It is a slightly different way to find a PDF compared to 
that of llliev et al.l (l2006al ) in which they sample independent sub-cubes in their simulation 
box. From these sampled spheres, we estimate mean (//), variance (<r 2 ), skewness (71), and 
kurtosis (72). We use three different sphere sizes of comoving radius 5, 10, and 20 Mpc/h, 
respectively. Figure 13 shows measurements of the underlying HII probability distribution 
functions. First of all, we see that the mean ionization fraction reaches 50% at z ~ 9.5 
(hereafter, z 50 %). Because the size of sampled spheres is much larger than the characteristic 
bubble sizes, the measured [i reflects the global mean ionization fraction. The largest dis- 
persion of ionization fraction appears at z ~ 8.9, close to z 50 %. It is easily seen that the 
degree of non-Gaussianity decreases with increasing size of the sphere, as shown in all three 
bottom panels, as expected. While the variance starts to fall at z < 8.9, the non-Gaussianity 
as measured by 71 and 72 continue to rise, peaking at z ~ 5.7 when the largest bubble starts 
to fill the entire simulation box. Interestingly, it might be of some comfort for analytic 
modelers that at z 50 %, the HII field is weakly non-Gaussian with 71 = (0.22,0.17,0.01) and 
72 = (—0.66, —0.26, —0.46), for r = (5, 10,20) Mpc/h, respectively. But 71 and 72 are most 
close to at z ~ 9 and ~ 10, respectively. Because the boundary between fully neutral 
and fully ionized regions is much more smaller than t he sizes of the sampling spheres, the 



meas ured non-Gaussianity is quite small around z 50 % (jWyithe &: Morales! 120071 : iLidz et al. 
2007h . 



3.4. Power spectrum of ionized hydrogen density fluctuation 

We present the 3D power spectrum of ionized hydrogen density fluctuation field (Snu), 
neutral hydrogen density fluctuation field (<5hii), and matter density fluctuation field (<5 mat ) 
in Figure 14. We note that the power spectrum is not measured for ionization fraction 
fluctuation but for ionized density fluctuation. The dimensionless power spectrum is defined 
here as 

= 7^y 47rk3 p ( k )- ( 6 ) 

We see that the ionized hydrogen density field shows a strong bias with respect to the 
matter power spectrum, where the bias is the strongest on small scales and decreases toward 
large scales. In general, the bias of ionized regions always decreases with time, whereas the 
bias of the neutral regions does the opposite. We find that Ajjj matches to that of 5 mat at 
the early time when most hydrogen is still neutral, as expected, whereas Aj^ n becomes equal 
to that of <5 ma t when the universe gets completely ionized. This is in agreement with the 
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underlying physics that the earlier bubbles are produced by highly biased galaxies formed 
inside DM halos of high a peaks, whereas only those highly biased hence large galaxies host 
most of the neutral hydrogen when the universe becomes highly ionized. It is interesting to 
note that at around z 50 % A^j and A^ n appear to have roughly the same shape, although the 
latter has somewhat more power than the former, while both have more power than the total 
matter, suggesting a significant non-gaussian distribution of both. We also note that there 
are significant difference between the nonlinear (actual) power spectrum of total matter and 
the linear power spectrum at k > 0.3 at z < 10, indicating the necessity to take nonlinear 
effects into account. 



4. CONCLUSIONS AND DISCUSSION 

Using the state-of-the-art, largest cosmological simulation of box size 100 Mpc/h with 
detailed radiative transfer of ionizing photons, we compute the evolution of HII bubbles 
during cosmological reionization from z ~ 25 to z ~ 6. Our simulation resolves galaxy 
sources that produce most of the ionizing photons in the universe in the concerned redshift 
range. Here are a few major findings from our analysis. 

(1) We find that, for significant HII bubbles, their number distribution is peaked at a 
volume of ~ 0.6 Mpc 3 /h 3 at all redshifts. But, at z < 10, one large, connected network of 
bubbles dominates the entire HII volume. 

(2) HII bubbles are highly non-spherical. This result is not totally unexpected, since 
one is quite familiar with the generic filamentary nature of cosmic structure formation. 
The mergers of HII bubbles blown by galaxies formed along filaments consequently produce 
filamentary (fatter) HII bubbles. 

(3) The HII regions are highly biased with respect to the underlying matter distribution 
with the bias decreasing with time as more and more less biased structures begin to dominate 
the ionizing photon production rate. The opposite is true for the neutral hydrogen region in 
the sense that the bias of the neutral regions increases with time. 

(4) The universe becomes 50% ionized at redshfit z ~ 9.5, when the HII region is 
actually the least non-gaussian in the redshift range concerned. The non-gaussianity of the 
HII region reaches its maximal near the end of the reionization epoch z ~ 6. But at all 
redshifts of interest there is a significant non-gaussianity in the HII field. 

(5) Population III galaxies play a very important role in the reionization process. With 
our spatially resolved treatment of metal enrichment, we show that small bubbles are initially 
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largely produced by Pop III stars. At z > 10 even the largest HII bubbles have a balanced 
ionizing photon contribution from Pop II and Pop III stars. At z < 8, however, Pop II stars 
start to dominate the overall ionizing photon budget for large bubbles, although Pop III 
stars continue to make a non- negligible contribution. 

(6) The relationship between halo number density and bubble size is complicated, al- 
though there is a strong correlation for large bubbles in which larger bubbles tend to have 
higher halo number density. In addition, we find that only the large bubbles (size > 1 
Mpc/h) contain halos of mass in excess of lO 9 M /h. 

Some of the results depend on modeling the formation of Pop III stars a nd the transition 



from Pop HI to Pop II star form ation that needs more careful studies (e.g. ISchneider et al. 



20061 ; ISmith fc Sigurdssonl 120071 ). First of all, the correct ionizing photon generation rate 
from the unit mass of Pop III stars or Pop II stars is as important as a careful treatment of 
recombination process in reionization simulation. As shown in §3.1.4, the formation of small 
bubbles is closely connected to the ionizing efficiency of the Pop III stars that is subject 
to the uncertain IMF of the Pop III stars. Our reionization simulation is the first trial to 
simulate the transition of Pop III star formation to Pop II star formation in a large simulation 
box of 100 Mpc/h. Im proving our understanding of chemical enrichment during reionization 



l^e.g. iGreif et al.l 120071 ). the transition of a dominant radiation source will be more certainly 
simulated in future. 

In addition to the uncertain physics of Pop III population, the numerical treatment 
of self-shielding and recombination has to be improved in future simulations. We used 
information of simulation particles for calculating the recombination rate that has better 
resolution than grid-based approach. But out approach also has limitation of resolving 
small-scale physics like self-shielding. That must be improved in the reionization simulation 
with a sub-grid model. 
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A. RADIATIVE TRANSFER: PHOTON- ADVECTION VERSUS 

RAY-TRACING 



We compare the photon-advection scheme with the ray-tracing scheme of iTrac &: Cen 



( 120061 ) by applying these RT algorithms to two reionization simulations with identical initial 
conditions and source prescriptions. The test is conducted in a 25 Mpc/h box with 720 3 dark 
matter particles and 90 3 RT cells and has the same spatial, mass, and temporal resolution 
as the large 100 Mpc/h simulation. 

We find that the photon-advection scheme produces very similar results in terms of the 
spatial and temporal evolution of HI for a substantial majority of the reionization epoch. 
Figure 15 shows very good agreement in the redshift evolution of the volume-averaged HI 
fraction /hi and deviations are only found when the box is already significantly ionized. 
The timing starts to differ at late stages near complete overlap. There are small delays of 
Az ~ 0.1 and 0.2 when the neutral fraction drops to 0.1 and 0.01, respectively. 

In order to quantify the agreement in the spatial evolution of HI, we cross- correlated 
the fm(x) fields from the two RT simulations and plotted the results in Figure 16. From the 
two power spectra, -P a dv(&) and P ray (fc), and the cross-power spectrum P a d v -ra y (&), we can 
define the bias function, 

" {k) s \SI (ai> 

and cross-correlation function, 

r(k) = JjgjXz£gxjj (A2) 
V^adv(fc)P ray (A;)' 

and use them to quantify the statistical correlation between the two fields. In both simula- 
tions, the neutral fraction drops to 0.5 and 0.25 at z — 8.2 and 7.5, respectively. At these two 
stages of reionization, the bias b(k) is very close to unity and only deviates by a maximum of 
~ 0.1 near the RT grid Nyquist frequency k = 11.3 h/Mpc. Similarly, the cross- correlation 
r(k) shows very good agreement even down to the smallest scale. In the ray-tracing sim- 
ulation, the neutral fraction drops to 0.1 at z = 7.1, but /hi is slightly higher by 0.04 in 
the photon-advection case. At this redshift, the bias and cross-correlation functions show 
that the two fields differ appreciably at all scales. However, when the two simulations are 
compared at the same neutral fraction of 0.1, the correlation is much better, particularly on 
large scales. The deviations from unity at the smallest scales are due to the appearance of 
new sources in the photon-advection simulation taken at a slightly later redshift. 

In summary, we find that the photon-advection scheme correctly captures the reion- 
ization process up until it is ~ 75% completed by volume. At earlier stages, the radiation 
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field is highly nonuniform, even within the HII regions, and the propagation of photons in 
the direction of decreasing radiation flux is a good description. However, at later stages of 
reionization near complete overlap, the radiation field is much more uniform and the weak 
gradients in the radiation density do not provide accurate directions for photon propaga- 
tion. We conclude that the photon-advection scheme provides a cost-effective approach to 
radiative transfer for a significant majority of the reionization epoch. However, the stages 
just before and after reionization should be simulated using more accurate approaches like 
ray-tracing. Our results in this paper are valid at all stages when the the ionization fraction, 
rather than the redshift, is used as an indicator of the progress of reionization. 
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Fig. 1. — Dark matter halo mass functions. Dark matter halos are identified using a friends- 
of-friends algorithm with a standard linking length of .2 times the mean inter-particle 
spaci ng;. Our mass func ti ons a re in better agreement with IWarren et al.l (120061 ) (solid) than 
with iPress fc Schechterl ( 119741 ) (dash) for z < 10. At z > 15, we under-resolve the halos 
because of the late starting redshift of z = 60. A smaller simulation starting at z = 300 
correctly captures the abundance of high redshift halos ( open circles ). 
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Fig. 2. — Comoving star formation rate (SFR). The total SFR (solid), from Pop III (long 
dash) and Pop II (dash) is consistent with that of iHernquist fc Springell (120031 ) (dotted). 
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Fig. 3. — Distributions of ionization fraction at redshifts z ~ 13.5, 11.2, 9.3, and 7.2. Highly 
ionized regions are represented by black while ionization fraction below 50 percent is shown 
as white. This is a plot of one slice in the simulation box. Each side of the plot is 100 Mpc/h. 
The global ionization fractions are ~ 10, 30, 50, and 90 percents for z ~ 13.5, 11.2, 9.3, and 
7.2, respectively. 
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Fig. 4. — Number and accumulated volume fraction distributions of HII bubble volumes 
(left) and bubble size distribution (right). At redshifts z = 20.6, 17.2, 13.5, 8.9, 6.3, and 5.7, 
respectively, the number (solid) and accumulated volume distribution (dash) is given with 
the total number of bubbles, N to t, and the global volume- weighted mean ionization fraction, 
< x >. Before large bubbles that are comparable to the simulation box size appear around 
the redshift 6, bubble size distributions show three peaks that have the volumes of about 0.6, 
0.03, and 0.006 Mpc 3 /h 3 . The domination of one large bubble appears at z < 10 as shown 
in the accumulated volume fraction s of bubbles . In t he right plot, we present the bubble 
size distribution that is measured as IZahn et al.l (120071 ) at z = 13.2, 12.1, 10.5, 9.0, and 8.0 
that are correspondent to the global volume-weighted mean ionization fraction of 0.12, 0.18, 
0.32, 0.52, and 0.75 from the small to large peak size. 
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Fig. 5. — Bubble merger history. The evolution of the bubble size is given for four positions 
that are occupied by the four randomly selected smallest bubbles at z = 21.5. The bubbles 
that occupy the four positions are tracked and their volumes are measured. This bubble 
merger history is well explained by the change of a dominant stellar population and its 
ionizing photon production rate between z ~ 10 and ~ 15 as shown in Figure 7. During that 
transition bubbles experience reconnection to other bubbles. 
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Fig. 6. — Evolution of ionizing photon production rate from Pop III and Pop II stars. The 
ionizing photon production rate is estimated by a star forming rate times an ionizing photon 
of hu > 13.6eV production rate per unit stellar mass. 
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Fig. 7. — Change of the number and volume fraction of HII bubbles, (upper) The total 
number of the HII bubbles increases until one large bubble dominates. Although the volumes 
of the smallest bubbles do not occupy a large fraction of the total volume, their number 
fraction dominates other sizes and shows a slightly later decline, (lower) Although the 
number of bubbles shows decrease after z ~ 10, the total volume fraction of all HII bubbles 
increases as cosmic reionization continues. The volume fraction of the largest bubble at 
each epoch explains the change of the bubble number is mainly caused by the fact that 
intermediate-size bubbles are merged to the largest bubble so that even the second largest 
bubble occupy a small fraction of the total volume. 
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Fig. 8. — Number density distribution of dark matter halos inside HII bubbles at z ~ 8 and 
10. Open triangle (red), dot (black), and open rectangle (blue) represent maximum, average, 
and minimum number density of halos for a given bubble size, respectively. The ranges of DM 
halo masses are M halo < 10 8 M /h, 10 8 M /h < M halo < 10 9 M /h, and 10 9 M /h < M halo . 
[See the electronic edition of the Journal for a color version of this figure.] 
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Fig. 9. — Ionizing photon production rate per unit volume inside HII bubbles at redshifts 
8.3 and 10.2. The ionizing photon production rate is calculated as used in Figure 6. Open 
triangle (red), dot (black), and open rectangle (blue) data points are maximum, average, and 
minimum values for different sizes of bubbles, respectively. Here we plot only bubbles that 
have stellar masses inside them. Note that ionizing photon production rate of Pop II stars 
is comparable to that of Pop III stars for bubbles of V > 0.6 Mpc 3 /h 3 . [See the electronic 
edition of the Journal for a color version of this figure.} 
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Fig. 10. — Distribution probability of the HII bubble shapes expressed as e\ and ei at z = 
13.8 and 6.3. Most bubbles are far from a sphere, i.e. e\ and €2 ^> 0. This trend does not 
change even though reionized volume becomes increasing. The distribution is derived from 
Equation 5. [See the electronic edition of the Journal for a color version of this figure] 
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Fig. 11. — Distribution probability of bubble shape parameters and volume at z = 13.8 and 
10. The bubble shape has a distribution that does not change even though the fraction of 
reionized regions increases between z = 13.8 and 10. While small bubbles have the narrow 
range of shapes around e% or e 2 ~ or 1, the bubbles of intermediate size show the broad 
range of shapes. The star symbols represent the shape of the largest bubble. [See the 
electronic edition of the Journal for a color version of this figure.} 
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Fig. 12. — Projected mean ionization fraction distribution of a single bubble at z = 13.8. 
Its derived shape parameters, e\ ~ 0.93 and €2 ~ 0.5, are consistent with this ionization 
fraction distribution. Here we set the background around the bubble to have zero ionization 
fraction for visualization. 
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Fig. 13. — Non-Gaussianity of reionization field. The statistics of mass- weighted ionization 
fraction are given for mean, variance, skewness, and kurtosis from top to bottom. Each 
statistic is calculated by sampling 5000 spheres of the comoving radius 5 (dot line, red), 10 
(dash line, green), and 20 (solid line, blue) Mpc/h, respectively. [See the electronic edition 
of the Journal for a color version of this figure.] 
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Fig. 14. — Power spectrum of 5m, 5un, and 5 ma t- At high redshift the power spectrum of a 
neutral hydrogen (dot line, red) is well matched to that of matter (dash line, green) while an 
ionized hydrogen (dot-dash line, blue) shows the match to the matter power spectrum at low 
redshift. A linearly developed matter power spectrum (solid line, black) does not match to 
a matter power spectrum at a small scale as a non-linear structure forms. [See the electronic 
edition of the Journal for a color version of this figure.] 



-35 - 




10 



15 



20 



Fig. 15. — A comparison of the volume-weighted neutral hydrogen fraction /hi from the 
photon-advection scheme (blue, solid) and ray-tracing scheme (red, dashed) for radiative 
transfer. The photon-advection scheme correctly captures the reionization process up until 
it is ~ 75% completed by volume. The magnitude of the difference in /hi (black, dotted) is 
generally very small and only reaches a maximum value of 0.05 near complete overlap. [See 
the electronic edition of the Journal for a color version of this figure] 
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Fig. 16. — A comparison of the photon-advection scheme and ray-tracing scheme for radiative 
transfer using the bias b(k) and cross-correlation r(k) of the neutral hydrogen fraction fm{x) 
field. At z — 8.2 (red, solid) and 7.5 (green, long-dashed), both simulations have volume- 
averaged neutral fractions of 0.5 and 0.25, respectively, and the spatial evolution of HI is 
very similar. At = 7.1 (blue, short-dashed), the correlation is poor because reionization has 
progressed slightly faster in the ray-tracing scheme. However, when the two simulations are 
compared at the same neutral fraction of 0.1 (magenta, dotted), the correlation is much 
better, particularly on large scales. [See the electronic edition of the Journal for a color 
version of this figure] 



